{
 "cells": [
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "import sisl\n",
    "import numpy as np\n",
    "import matplotlib.pyplot as plt\n",
    "from functools import partial\n",
    "%matplotlib inline"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "TBtrans is capable of calculating transport in $N\\ge 1$ electrode systems. In this example, we will explore a 4-terminal graphene nano-ribbon (GNR) crossbar (one zig-zag GNR, the other armchair GNR) system."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "graphene = sisl.geom.graphene(orthogonal=True)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "R = [0.1, 1.43]\n",
    "hop = [0., -2.7]"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Create the two electrodes in $x$ and $y$ directions. We will force the systems to be nano-ribbons, i.e. only periodic along the ribbon. In `sisl` there are two ways of accomplishing this.\n",
    "\n",
    "1. Explicitly set number of auxiliary supercells\n",
    "2. Add vacuum beyond the orbital interaction ranges\n",
    "\n",
    "The below code uses the first method. \n",
    "\n",
    "Please see if you can change the creation of `elec_x` by adding vacuum.  \n",
    "**HINT**: Look at the documentation for the `sisl.Geometry` and search for vacuum. To know the orbital distance, look up `maxR` in the geometry class as well."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "elec_y = graphene.tile(3, axis=0)\n",
    "elec_y.set_nsc([1, 3, 1])\n",
    "elec_y.write('elec_y.xyz')\n",
    "elec_x = graphene.tile(5, axis=1)\n",
    "elec_x.set_nsc([3, 1, 1])\n",
    "elec_x.write('elec_x.xyz')"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Subsequently, we create the electronic structure."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "H_y = sisl.Hamiltonian(elec_y)\n",
    "H_y.construct((R, hop))\n",
    "H_y.write('ELEC_Y.nc')    \n",
    "H_x = sisl.Hamiltonian(elec_x)\n",
    "H_x.construct((R, hop))\n",
    "H_x.write('ELEC_X.nc')"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Now we have created the electronic structure for the electrodes. All that is needed is the electronic structure of the device region, i.e. the crossing nano-ribbons."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "dev_y = elec_y.tile(30, axis=1)\n",
    "dev_y = dev_y.translate( -dev_y.center(what='xyz') )\n",
    "dev_x = elec_x.tile(18, axis=0)\n",
    "dev_x = dev_x.translate( -dev_x.center(what='xyz') )"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Remove any atoms that are *duplicated*, i.e. when we overlay these two geometries, some atoms are the same."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "device = dev_y.add(dev_x)\n",
    "device.set_nsc([1,1,1])\n",
    "duplicates = []\n",
    "for ia in dev_y:\n",
    "    idx = device.close(ia, 0.1)\n",
    "    if len(idx) > 1:\n",
    "        duplicates.append(idx[1])\n",
    "device = device.remove(duplicates)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Can you explain why `set_nsc([1, 1, 1])` is called? And if so, is it necessary to do this step?\n",
    "\n",
    "---\n",
    "\n",
    "Ensure the lattice vectors are big enough for plotting.\n",
    "Try and convince yourself that the lattice vectors are unimportant for TBtrans in this example.  \n",
    "*HINT*: what is the periodicity?"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "device = device.add_vacuum(70, 0).add_vacuum(20, 1)\n",
    "device = device.translate( device.center(what='cell') - device.center(what='xyz') )\n",
    "device.write('device.xyz')"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Since this system has four electrodes, we need to tell TBtrans where the 4 electrodes are in the device. The following lines prints out the fdf-lines that are appropriate for each of the electrodes (`RUN.fdf` is already filled correctly):"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "print('elec-Y-1:   semi-inf -A2: {}'.format(1))\n",
    "print('elec-Y-2:   semi-inf +A2: end {}'.format(len(dev_y)))\n",
    "print('elec-X-1:   semi-inf -A1: {}'.format(len(dev_y) + 1))\n",
    "print('elec-X-2:   semi-inf +A1: end {}'.format(-1))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "H = sisl.Hamiltonian(device)\n",
    "H.construct([R, hop])\n",
    "H.write('DEVICE.nc')"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Exercises\n",
    "\n",
    "In this example, we have more than one transmission path. Before you run the below code which plots all relevant transmissions ($T_{ij}$ for $j>i$), consider if there are any symmetries, and if so, determine how many different transmission spectra you should expect? Please plot the geometry using your favourite geometry viewer (`molden`, `Jmol`, `gdis`, ...). The answer is not so trivial."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "tbt = sisl.get_sile('siesta.TBT.nc')"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Make easy function calls for plotting energy resolved quantities:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "E = tbt.E\n",
    "Eplot = partial(plt.plot, E)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Make a shorthand version for the function (simplifies the below line)\n",
    "T = tbt.transmission\n",
    "t12, t13, t14, t23, t24, t34 = T(0, 1), T(0, 2), T(0, 3), T(1, 2), T(1, 3), T(2, 3)\n",
    "Eplot(t12, label=r'$T_{12}$'); Eplot(t13, label=r'$T_{13}$'); Eplot(t14, label=r'$T_{14}$');\n",
    "Eplot(t23, label=r'$T_{23}$'); Eplot(t24, label=r'$T_{24}$'); \n",
    "Eplot(t34, label=r'$T_{34}$');\n",
    "plt.ylabel('Transmission'); plt.xlabel('Energy [eV]'); plt.ylim([0, None]); plt.legend();"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "- In `RUN.fdf` we have added the flag `TBT.T.All` which tells TBtrans to calculate *all* transmissions, i.e. between all $i\\to j$ for all $i,j \\in \\{1,2,3,4\\}$. This flag is by default `False`, why?\n",
    "- Create three plots each with $T_{1j}$ and $T_{j1}$ for all $j\\neq1$."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Insert plot of T12 and T21"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Insert plot of T13 and T31"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Insert plot of T14 and T41"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "- Considering symmetries, try to figure out which transmissions ($T_{ij}$) are unique?\n",
    "- Plot the bulk DOS for the two differing electrodes.\n",
    "- Plot the spectral DOS injected by all four electrodes."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Helper routines, this makes BDOS(...) == tbt.BDOS(..., norm='atom')\n",
    "BDOS = partial(tbt.BDOS, norm='atom')\n",
    "ADOS = partial(tbt.ADOS, norm='atom')"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Bulk density-of-states:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "Eplot(..., label=r'$BDOS_1$');\n",
    "Eplot(..., label=r'$BDOS_2$');\n",
    "plt.ylabel('DOS [1/eV/N]'); plt.xlabel('Energy [eV]'); plt.ylim([0, None]); plt.legend();"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Spectral density-of-states for all electrodes:\n",
    "- As a final exercise, you can explore the details of the density-of-states for single atoms. Take for instance atom 205 (204 in Python index) which is in *both* GNRs at the crossing.  \n",
    "Feel free to play around with different atoms, subset of atoms (pass a `list`), etc.  "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "Eplot(..., label=r'$ADOS_1$');\n",
    "...\n",
    "plt.ylabel('DOS/N [1/eV]'); plt.xlabel('Energy [eV]'); plt.ylim([0, None]); plt.legend();"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "- For 2D structures one can easily plot the DOS per atom via a scatter plot in `matplotlib`. Here is the skeleton code for that. You should select an energy point and figure out how to extract the atom resolved DOS (you will need to look up the documentation for the `ADOS` method to figure out which flag to use)."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "Eidx = tbt.Eindex(...)\n",
    "ADOS = [tbt.ADOS(i, ....) for i in range(4)]\n",
    "f, axs = plt.subplots(2, 2, figsize=(10, 10))\n",
    "a_xy = tbt.geometry.xyz[tbt.a_dev, :2]\n",
    "for i in range(4):\n",
    "    A = ADOS[i]\n",
    "    A *= 100 / A.max() # normalize to maximum 100 (simply for plotting)\n",
    "    axs[i // 2][i % 2].scatter(a_xy[:, 0], a_xy[:, 1], A, c=\"bgrk\"[i], alpha=.5);\n",
    "plt.xlabel('x [Ang]'); plt.ylabel('y [Ang]'); plt.axis('equal');"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python 3",
   "language": "python",
   "name": "python3"
  },
  "language_info": {
   "codemirror_mode": {
    "name": "ipython",
    "version": 3
   },
   "file_extension": ".py",
   "mimetype": "text/x-python",
   "name": "python",
   "nbconvert_exporter": "python",
   "pygments_lexer": "ipython3",
   "version": "3.11.6"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 4
}
